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Abstract 

Inhomogeneity, in its many forms, appears frequently in practical physical systems. Readily apparent in quantum 
systems, inhomogeneity is caused by hardware imperfections, measurement inaccuracies, and environmental variations, 
and subsequently limits the performance and efficiency achievable in current experiments. In this paper, we provide 
a systematic methodology to mathematically characterize and optimally manipulate inhomogeneous ensembles with 
concepts taken from ensemble control. In particular, we develop a computational method to solve practical quantum 
pulse design problems cast as optimal ensemble control problems, based on multidimensional pseudospectral approx- 
imations. We motivate the utility of this method by designing pulses for both standard and novel applications. We 
also show the convergence of the pseudospectral method for optimal ensemble control. The concepts developed here 
are applicable beyond quantum control, such as to neuron systems, and furthermore to systems with by parameter 
uncertainty, which pervade all areas of science and engineering. 

I. Introduction 

Recent advancements in quantum research have enabled breakthroughs in biology, chemistry, physics, engineering, 
and medicine including better methods to understand the structure of macromolecules used in biochemical signaling 
and drug delivery, to facilitate the fast and efficient storage of information, and to yield higher resolution medical 
images for diagnosis and treatment of early stage cancer |IT]-|l3|. Most, if not all, measurements and manipulations 
of quantum systems are achieved through the appropriate design of externally apphed time-varying electromagnetic 
pulses, or controls |[T]. These pulses guide the system to produce a desired time-evolution or a specific terminal 
state. The design of such pulses is made significantly more difficult by inherent variations within the systems of 
interest. Inhomogeneity is one of the fundamental obstacles for the practical implementation and physical realization 
of quantum science and quantum technology. In classical systems the dispersions resulting from inhomogeneity is 
often compensated for by feedback control. Significant research effort has been employed in the area of quantum 
feedback control with several promising theoretical and practical discoveries in recent years |@], lH]. There is still 
a large portion of quantum systems for which state feedback is either impractical or difficult to achieve due to 
the short timescales and large state-space of quantum phenomena. These limitations motivate us to consider the 
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open-loop synthesis of optimal pulses that achieve a desired goal while compensating for the inhomogeneity present 
in the quantum ensemble. 

The behavior of a bulk quantum system is the aggregate behavior of a large ensemble of individual quantum 
systems. Although in isolation these individual systems, e.g. atoms, spins, qubits, etc., are fundamentally identical, in 
a physical system they are distinct due to different chemical and electromagnetic environments. This variation across 
the ensemble exhibits itself at the macroscopic level as variation in the values of parameters that characterize the 
dynamics of the bulk quantum system |6|. For example, adjacent atoms within the same macromolecule shield the 
full strength of the applied external magnetic fields. Varied levels of shielding create a dispersion in the frequency 
of the quantum spins, which is observed as inhomogeneity in the value of the natural frequency of the bulk sample. 
In addition, hardware imperfection causes attenuation in the applied electromagnetic pulse over the ensemble and 
can be represented as variation in a scaling factor multiplying the applied pulse. Often several pulses are applied 
in sequence in order to achieve an intricate time-evolution of the system Q. Each pulse is designed assuming an 
exact (usually uniform) initial system state, however, in practice, the previous pulses only prepare the system to 
within a neighborhood of the assumed initial state. The additive error in such a pulse sequence can cause significant 
performance degradation. 

Guiding the evolution of inhomogeneous ensembles is a central idea in the design and implementation of quantum 
experiments. As such, there is a rich literature of methods addressing this class of challenging problems. Initially 
these were intuitive or ad-hoc methods motivated by the symmetry of the state space Q, ISj, which were then 
augmented with various heuristic and specifically designed techniques ||9l, lITOl . More recently pulse design problems 
have been cast as optimal control problems ifTTl - lfTSl . Here we present a methodology that addresses the difficulties 
of the current methods and is easily generalizable to any inhomogeneous ensemble or uncertain system. The proposed 
method has both theoretical, such as convergence rates and computational complexity, and practical, such as ease 
of implementation and computation time, advantages. 

In this article we describe a framework to pose robust quantum pulse design problems in the language of mathe- 
matical control theory with support from new theoretical concepts in ensemble control lfT6l - lfT8l and computational 
advances in multidimensional pseudospectral methods adapted for ensemble systems ||6l, lfT9l . In a larger context, 
we provide a rigorous methodology to study and control inhomogeneous ensembles or systems with parameter 
uncertainty from any area or application. In the following section we introduce the problem statement as well as 
our theoretical and computational tools. In Section Hill we take several examples from nuclear magnetic resonance 
(NMR) in liquids modeled by the bilinear Bloch equations, including broadband excitation in the presence of 
inhomogeneity, a sequence of broadband pulses robust to variation in the initial conditions, and systems with a 
time-varying frequency. We then provide empirical and theoretical justifications that the solutions computed using 
the pseudospectral method converge to solution of the original optimal control problem. 
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II. Motivation & Theory 

In this section, we review the underlying concepts involved in our approach to design robust quantum pulses 
as well as the broader mathematical formulations necessary to characterize and solve such design problems. In 
what follows, we present a highly general model of quantum dynamics, which demonstrates the abundance of 
inhomogeneity in these problems and motivates studying the control of a family of parameterized systems. We then 
show how the notion of ensemble control is well suited for dealing with the inherent variation and uncertainty in 
practical quantum systems and formulate a new type of optimal control problem based on ensemble control. 

A. Quantum Dynamics & Pulse Design 

The dynamics of a quantum system is given by the time-evolution of its density matrix. We consider here general 
dynamics in which the system may have interaction with the environment that leads to dissipation in the system 
state. Under the Markovian approximation, where the environment is modeled as an infinite thermostat which has 
constant state, the evolution of the density matrix can be written in Lindblad form in terms of the system Hamiltonian 
H{t) and superoperator L{-) which model the unitary and nonunitary dynamics [20 J . respectively. 

The expression of the Hamiltonian has components corresponding to free evolution Hamiltonian, Hj, and the control 
Hamiltonians Hi, 

m 

i=l 

where Ui{t) are externally applied electromagnetic pulses that can be used to manipulate, or guide, the evolution of 
the system state. Typical pulse design problems involve designing these pulses, or controls, to bring the final state 
of the density matrix p{T) as close as possible to a target operator This problem can be transformed, by taking 
the expectation values of the operators involved in the state transfer, to the vector-valued, bilinear control problem, 

X e R" and u e M™ given by. 



d 
dt 



(1) 



where Hd £ M"^" corresponds to the drift evolution representing Hf and L, Hi G R"^" corresponds to the 
controlled evolution representing Hi, and t G [0, T] ||2T1 . While ([T]i accurately represents the classical interaction of 
magnetic fields, in practice the effective fields - and, therefore, the matrices representing the Hamiltonians Hd and 
Hi - show variation in magnitude due to different chemical environments and equipment errors. The system can no 
longer be described by a single equation but rather by a family of equations with variation in the parameters that 
characterize the motion, which motivates us to consider the dispersion in the dynamics as a continuum parameterized 
by the system values, 

d 

-x{t, s) ^ [Hd{s) + u,{t)H^{s)\ x{t, s), (2) 
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where s G C M"* is a d-dimensional interval representing the d parameters exhibiting variation |[T9]- In a 
more general formulation the matrices representing the Hamiltonians can be time-dependent, Hd = s) and 

Hi = T~ii{t, s), as in the case of random fluctuations. 

Designing a single set of controls (pulses) Ui{t) that simultaneously steer an ensemble of dispersive systems 
in (|2|i from an initial state to a desired final state is a fundamental problem in the control of quantum systems. 
Moreover, similar parameterized structures can be found across all areas of science and engineering, such as in 
neuroscience where a single stimulus is used to trigger a simultaneous firing of neuron oscillators with distinct 
oscillation frequencies ll22l . In these applications full state feedback, which is required in most current methods to 
compensate for system uncertainty, is impractical to obtain due to the sheer number of members (and states) within 
the ensemble. Averaged measurement is possible in some applications, however, this type of measurement restricts 
the forms of available feedback. It is, then, of particular importance to consider the corresponding open-loop control 
problem. 

B. Optimal Ensemble Control 

Systems as in dU motivate the study of a new class of inhomogeneous control systems. Ensemble control ifTTl 
is a mathematical framework to characterize parameterized systems of the form. 



where x G R", u G M"*, s G -D C M'^, with F and xo{s) smooth functions of their respective arguments. The 
significant challenge of this class of control problems originates from requiring the same open-loop control, u{t) 
to guide the continuum of systems from an initial distribution, xq{s), to a desired final distribution, over the 
corresponding function space. Fundamental properties of these systems, such as controllability, are of particular 
interest - specifically addressing what types of inhomogeneities can be compensated for robustly. For example, it 
has been shown that the controllability of an ensemble of bilinear Bloch equations, used as a sample system in this 
paper (see Section Ulll i. corresponds to the synthesis of appropriate polynomials lfT6l and controllability conditions 
for an ensemble of time-varying linear systems are related to the Picard criterion of Fredholm integral equations 
of the first kind (TT\. 

Subsequently, given dynamics and initial and final distributions, we seek methods to construct controls for such 
steering problems. As with any control problem, in general there may be many possible solutions that satisfy a state- 
to-state ensemble control problem. In addition there are often benefits, penalties, and limitations that are associated 
with the physical system, which can be used to rank the different solutions. Such practical considerations lead to 
considering an optimal control problem based on the ensemble dynamics in (O which includes a cost functional 
(with terminal, (p, and running, C, cost terms) to be minimized as well as possible endpoint and path constraints 





(3) 
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(e and g, respectively), 

ip{T,x{T,s))+ [ C{x{t,s),u{t))dt ds, (4) 

D Jo 

s.t. ^^(^1 s) — F(t, s, x{t, s), u{t)) , 
e{x{0,s),x{T, s)) = 0, 
g(a;(t,s),u(t)) < 0. 

An optimal nonlinear control problem of this form is, in general, analytically intractable. Computational methods are 
then required to solve such exceedingly complex optimal ensemble control problems. The idea from our previous 
work that constructing appropriate polynomials is a key tool in characterizing the controllability of ensemble 
systems of interest motivates the use of polynomials within the computational framework lfT6i . Below we review 
the main ideas of the previously established pseudospectral method for optimal control to lay the foundation for 
our developed extension to optimal ensemble control problems. In Section |IV] we complete this framework with a 
proof of convergence of this numerical method. 

Without loss of generality, we consider a general continuous-time optimal control problem defined on the time 
interval SI = [—1, 1], which can be achieved by a simple affine transformation. 

Problem 1 (Continuous -Time Optimal Control): 

min J{x,u) ^ ip{x{l)) + j C{x{t),u{t))dt, (5) 

s.t. ^xit) = fit,xit),u{t)), (6) 
dt 

e(x(-l),a;(l)) = 0, (7) 

gix{t),uit)) <0, (8) 

\Ht)\\^<A, ueH^in),a>2 (9) 

where G C" is the terminal cost; the running cost, C e C", where C" is the space of continuous functions with 
a classical derivatives, and dynamics, / G C"^^, where C"^^ is the space of n-vector valued C"~^ functions, 
with respect to the state, x{t) £ R", and control, u{t) E W"; e and g are terminal and path constraints, respectively; 
H"^{il) is the TO-vector valued Sobolev space. The norm associated with the Sobolev space with m = 1, H"{n), 
is given with respect to the L^(Sl) norm ||231 . 

^ k=0 

C. Pseudospectral Method 

The pseudospectral method was originally developed to solve problems in fluid dynamics and since then has been 
successfully used for optimal control Il24l - ll26l and applied to various areas 11211 . ||221 . Pseudospectral discretization 
methods use expansions of orthogonal polynomials to approximate the states of the system and thereby inherit the 
spectral accuracy characteristic of orthogonal polynomial expansions (the /c* coefficient of the expansion decreases 
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faster than any inverse power of k) [23|. Through special properties, derivatives of these orthogonal polynomials can 
be expressed in terms of the polynomials themselves, making it possible to accurately approximate the differential 
equation that describes the dynamics with an algebraic relation imposed at a small number of discretization points. 
An appropriate choice of these discretization points, or nodes, facilitates the approximation of the states as well as 
ensuring accurate numerical integration through Gaussian quadrature. 

As a collocation (or interpolation) method, the pseudospectral method uses Lagrange polynomials to approximate 
the states and controls of the optimal control problem, 

xit) « /jvx(t) = ^f^o Xfe4(t), (10) 

Uit) « iNuit) = ELo^kikit), (11) 

where x(tk) = lNx{tk) — Xk and u{tk) = lNu{tk) = Uk because the Lagrange polynomials have the property 
(k{ti) — Ski, where Ski is the Kronecker delta function and tk are the interpolation nodes ll27l . Therefore, the 
coefficients Xk and Uk are the discretized values of the original problem and become the decision variables of the 
subsequent discrete problem. 

Although the interpolation with Lagrange polynomials discretizes the original problem, we require a means to 
ensure that both the integral in the cost functional is computed accurately and the dynamics are obeyed. The integral 
can be approximated through Gauss quadrature; here we use Legendre polynomials as the orthogonal basis for the 
pseudospectral method. The Legendre-Gauss-Lobatto (LGL) quadrature approximation, 

£,{t)dt, (12) 



^1 1 

J ^ /(t)di « ^ /(t,)u;„ w,^ y ^. 



is exact if the integrand / G P2Ar-i and the nodes ti E p^'-^^, where P2Ar-i denotes the set of polynomials of 
degree at most 2N - 1 and where T^*^^ = {U : = 0, z = 1, . . . , iV - 1}U{-1. 1} the iV + 1 LGL 

nodes determined by the derivative of the iV* order Legendre polynomial, LN{t), and the interval endpoints ||231 . 

Using the LGL nodes, we can rewrite the Lagrange polynomials in terms of the orthogonal Legendre polynomials, 
which is critical to inherit the special derivative and spectral accuracy properties of the orthogonal polynomials 
despite using Lagrange interpolating polynomials. Given tk € F^'-^^, we can express the Lagrange polynomials as 

ED, 

e (t) - 1 jt' - i)LN{t) 

''^ ' N{N + l)LN{tk) t-tk ■ 
The derivative of ^ at tj e T^^^ is then, 

N N 

dt 

where D is the constant differentiation matrix 

We are now able to write the discretized optimal control problem using equations (fTol i. (fTTT i. (fT2l i. and (fTsl i. We 



^-jNx{tj) = ^ Xkikitj) = ^ DjkXk = {'DNx){tj), (13) 

k=0 k=0 
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transform the continuous-time problem to a constrained optimization, 

N 

min (p{xn) + ^^C{xi,Ui)wf' , 



i=0 

N 



s.t. ^DjkXk ^ f(xj,Uj), (14) 

fc=0 

e{xQ,XN) = 0, 

g{xj,u,)<0, Vj G {0,1,...,7V}. 

D. Multidimensional Pseudospectral Method 

The pseudospectral method lends itself to a natural extension to consider the ensemble case, which we develop 
here. This is most readily apparent for a single parameter variation, i.e., s E [a, b] C M, however, is easily scaled 
to an arbitrary parameter dimension. In this basic case, the ensemble extension of ( fTOl ) is 

N N / N^, \ 

X{t,s) lNxN,x{t,s) = '^Xk{s)ik{t) «^ '^XkrUis) ikif). (15) 



fe=0 fe=0 \r=0 



The approximate derivative from (ITJt at the LGL nodes in the respective t and s domains, G F^'-^^ and Sj G V^'^^ 



IS 



^lNy.N,x{ti,Sj) = ^ Afc ( ^a;fcr4(sj) I = ^ Afca^fcj, (16) 

A;=0 \r=0 / A;=0 

where Xfcj = a;(tfe, Sj). In these equations we use a two-dimensional interpolating grid at the + 1 and A^^ + 1 
LGL nodes in time and the parameter, respectively. For a general number of parameters, s = (si, S2, • ■ • , Sd)' G 
D C M'^,d > 1, 

AT N N,^ N,^ 

X{t,s) « lNxN,-^x---xN,^x{t,s) = ^ife(s)4(i) = X! X! ' ■ ' X! ^kri...rJrA^d) ' ■ ■ tri{si)tk{t) ■ (17) 

k=0 fc=Ori=0 rd=0 



and the derivative is, correspondingly, with j = (ji, j2, • ■ • , jd)', 

d ^ 

-^lNxN^-^x---xN^^xit, Sj) — '^^DikXkj^...j^- (18) 

fc=0 

The simplification from (fTTI i to (fTSl l illustrates why the pseudospectral approximations are effective methods for 
ensemble control, as they mimic the lack of information in the parameter dimension. This aspect will also make 
the extension of the convergence proof for ensemble systems straightforward, as will be discussed in Section |IV] 

III. Examples 

In this paper, we consider several examples based on the prototypical quantum control system described by 
the Bloch equations |30|. The Bloch equations have been found to model a range of quantum phenomena from 
protein spectroscopy in nuclear magnetic resonance (NMR) (|T| and medical scans in magnetic resonance imaging 
(MRI) [311 to Rabbi oscillations in quantum optics f32l. In the following discussion, we will consider the specific 
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application and terminology for NMR spectroscopy, however, the methods and results are easily transferred to these 
other areas of interest. In NMR spectroscopy, when the duration of the pulse design problem is small compared 
with the relaxation times (T ^ Ti, T2, the characteristic longitudinal and transverse relaxation times, respectively), 
the evolution of spins can be well approximated as sequences of unitary rotations driven by the static magnetic field 
and the applied electromagnetic controls. In practice, the effective fields generating these rotations show variation 
across the quantum sample due to hardware imperfection and chemical shielding, which leads us to consider a range 
of magnetic field variations. The corresponding dimensionless Bloch equations in the rotating frame (see Appendix 
|A]i are, 

^M{t,uj,e)= u}n,_ + eu{t)ny + ev{t)n^ M{t,uj,e), (19) 

where Af(t,aj,e) = (Mx{t,u!,e),My{t,u!,e),Mz{t,uj,e)) is the Cartesian magnetization vector for the parameter 
values s = (w,e), uj G [—B,B] C K, is the dispersion of natural frequencies, e G [1 — 6,1 + 6], < S < 1, 
is the amplitude attenuation factor, and fia G SO(3) is the generator of rotation around the a axis. A pulse that 
compensates for the dispersion in frequency and is insensitive to the scaling of the applied controls is called 
a broadband pulse robust to inhomogeneity. In this section we consider several examples based on this model, 
including pulses robust not only to frequency dispersion and inhomogeneity, but also robust to uncertainty in initial 
conditions and time-varying frequencies. 

A. Robust TT Pulse 

Variation and dispersion in system dynamics pervade all physical experiments. In quantum systems, these 
inhomogeneities are often large enough to cause significant reduction in performance. The systematic framework 
we present here provides a rigorous way to frame any general pulse design problem for quantum control, as well 
as other areas of parameterized and uncertain systems. 

A canonical problem in the control of quantum systems modeled by the Bloch equations is to design pulses that 
will accomplish a state-to-state transfer of the system. Such pulses, e.g., 7r/2 and tt pulses (accomplishing 7r/2 
and TT rotations, respectively), are the fundamental building blocks of the pulse sequences used in many quantum 
experiments. Here, consider the inversion, or tt, pulse that rotates the net magnetization from the equilibrium position 
(+z) to the —z axis, i.e., M(0) = (0 1)' M{T) ==(0 0— 1)'. In the ensemble case, this goal corresponds to 
a uniform inversion of the spin vector across all choices of frequency and inhomogeneity. Specifically we consider 
the optimal ensemble control problem. 



min / / M^{T,uj,e) dw de+ u'^{t) + v'^{t) dt, 

Jl-S J-B Jo 

s.t. ^Af(t,w,e)= ujQ^ + eu{t)ny + ev{t)i}^ M{t,uj,€), 
M(0,w,e) = (0 1)', 

Vu^{t)+v^{t) <A, yte [o,r], 



(20) 
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where A is the maximum allowable amplitude and the cost functional serves to equally minimize the z-component 
of the spin vector (integrated across the ensemble) and the energy of the designed pulse. 

Figure [T] displays a pulse that compensates for B = 1 and 5 = 0.1 (10%) as well as the corresponding inversion 
profile. In physical units for a normalizing amplitude of 10 kHz, the maximum amplitude is A = 20 kHz with 
bandwidth uj G [—20, 20] kHz and duration T — 120 /is. Pulses developed in this manner have been implemented 
experimentally in true protein NMR experiments to yield significant improvement in signal recovery [|6]. Although 
designing individual pulses is of importance and benefit there are a myriad of other variations and uncertainties 
within typical quantum experiments, which calls for an approach that can address such new inhomogeneities and 
their corresponding challenges. 

B. Uncertainty in Initial Conditions 

In most experiments, individual pulses, such as the one in Figure [T] are combined into a longer pulse sequence, 
which performs a more complicated manipulation of the system state with intermediate steps and goals. Even in 
the case of highly optimized individual pulses, as shown in the prior example, there is an error between the desired 
and actual final states. Moreover, pulses depend upon an exact (and usually uniform) initial condition in order to 
achieve their expected levels of performance. These effects combine to create a magnified accumulated error at the 
termination of the pulse sequence. The variation of the initial conditions of these pulses, therefore, causes significant 
degradation in achievable performance. 

A representative example of such a pulse sequence is to perform a three step pulse sequence, which rotates the 
magnetization of the ensemble (1) from equilibrium i+z) to a point on the transverse plane (e.g. +j/); (2) to the 
opposite point on the transverse plane (e.g. — y); (3) back to the equilibrium position (+z). Such pulses generally 
include "phase locking" pulses before and after the second pulse during which the magnetization dissipates. This 
dissipation is the portion of the experiment that is important to recover accurately and reflects a quantity to be 
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Fig. 2. Pulses are optimized to produce a desired z y ^ —y — > z evolution of the Bloch equations. The upper plot displays the concatenation 
of individually optimized z y and y — > —y pulses, which achieves the dashed terminal profiles shown below, with respective average 
performances: 0.99, 0.98, 0.97 (0.91 minimum). The middle plot displays a 3-part simultaneously-optimized pulse robust to variation in the 
initial condition and achieves the solid terminal profiles shown below, with respective average performances: 0.99, 0.99, 0.99 (0.97 minimum). 
The noticeable enhancement in perfonnance and uniformity is due to compensating for the inhomogeneity in the initial condition of the individual 
pulses. 



measured, for example, a metabolic rate ll33l . Il34l . If, in addition, there is accumulated error due to uncertainty in 
the initial conditions of the individual pulses, this leads directly to measurement inaccuracy. Here, by removing the 
"phase locking" pulses, we can abstract this pulse sequence to a unitary process and directly address any losses 
due to error The controllability of the Bloch equations is shown constructing parameter-dependent (e.g. frequency, 
rf inhomogeneity) rotations of the spin vectors [IS]. This, therefore, ensures that the problem with variation in 
initial conditions can be solved provided that the initial conditions can be parameterized by the frequency and 
rf-inhomogeneity. 

Figure 12] displays a three-stage optimized pulse designed by the multidimensional pseudospectral method which is 
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t t t 

Fig. 3. Control pulses (top) and state trajectories (bottom) coiTesponding to different objectives and designed to compensate for the time- 
varying frequency Lo{t) = sin(t). A single-system state transfer Af(0) = (0 1)' — > M(T) = (10 0)' is designed using the terminal cost 
^(T) = Mx{T) and running costs C{t) = (left), C(t) = 0.l{u{t)'^ + f (t)^) (middle), C(t) = 0.1 (right). The terminal time was free in 
all cases, bounded by Tmax = 1. 

robust to frequency dispersion and variation in the initial conditions of the three stages. This pulse was run as three 
concurrent optimizations, with the final states of one pulse fed in as the initial conditions of the next. This optimized 
pulse is compared with the combination of three separately optimized pulses; these combined pulses were designed 
with equal total duration. The terminal profiles at each intermediate goal quickly show the evidence of accumulated 
error in the case of the individually optimized pulses (each individual pulse has an average performance greater 
than 0.98). Most importantly, the uniformity of the inversion is lost in the additive error, with dips in performance 
down to 0.91. 

C. Time-Varying Frequency 

Until now, we have considered that the dispersion and uncertainty of the system are stationary. However, 
addressing time-varying fluctuations in parameters is also of particular theoretical and practical importance. For 
example, in the formulation of quantum control problems given in dU we noted that the Hamiltonians can be time- 
varying, motivated by such phenomena as random telegraph noise 1351 . The first step to addressing stochastic 
variations in such physical systems is to demonstrate control of time-varying systems, such as given by the 
expectation value of the corresponding random process. 

Figure |3] presents a series of optimizations designing 7r/2 pulses providing a state transfer +z to +x, while 
compensating for a time-varying frequency, Lu{t) — sm{t). Various choices of cost functional yield different results. 
The arbitrary control pulse profile corresponding to the terminal cost <p(T) = (T) (Fig. |3] left) motivates studying 
optimal control methods that provide the capacity for hybrid objectives resulting in more physically meaningful 
controls, e.g. minimizing energy (middle) and time (right). 
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IV. Convergence 

By accepting and implementing a numerical method, we implicitly assume that the transformations and discretiza- 
tion used to prepare the problem for computational work does not fundamentally alter the nature of the problem. It 
is then critically important to show that this assumption is justified. Here we do so by both empirical and theoretical 
means. More specifically, we show that as the number of discretizations in the pseudospectral method (and samples 
in the multidimensional pseudospectral method) increases the solution of the algebraic nonhnear programming 
problem converges to the solution of the original continuous-time optimal control problem. For this argument, we 
consider a modified nonlinear programming problem statement. 

Problem 2 (Algebraic Nonlinear Programming): 

N 

min J{x,u) = ip[xn) + '^C[xk,Uk)wk (21) 

fe=0 

s.t. \\f{XNX,XNu)~VNx\\j^<CdN^-°' (22) 

e(So,XAr) = (23) 

g{xk.uk)<0 (24) 

\\uk\\<A Vfc = 0,l,...,A^ (25) 

where is a positive constant; we define the discrete i^j(il) norm \\h\\is[ ~ ^ {h, ft.) at, for h,hi,h2 £ £^^(17), 
n = [-1, 1], with, 

N 

{hi,h2)N = y^^h[{tk)h2{tk)wk, 

where ' denotes the transpose and Wk is the Gauss quadrature weight from (fT2t . 

Remark 1: The dynamics in (l22T i have been relaxed from the equality in (fT4l i to ensure the feasibility of the 
discrete problem, which is used in Proposition [T] It is trivial to show that in the limit, as — s> oo, these two 
conditions coincide. 

We seek to address three questions related to solving the continuous-time optimal control (Problem [T) by solving 
the pseudospectral discretized constrained optimization (Problem |2]i. Suppose a feasible solution {x, u) exists to 
Problem [T] Under what conditions: 

1) Feasibility: For a given order of approximation, N, does Problem |2] have a feasible solution, [x, u), which are 
the interpolation coefficients given in (fTOl i and (fTTli? 

2) Convergence: As N increases, does the sequence of optimal solutions, {{x\u'')}, to Problem [2] have a 
corresponding sequence of interpolating polynomials which converges to a feasible solution of Problem \Q 
Namely, 

lim {Tnx\Inu^) — {x,u) 

N~foo 
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3) Consistency: As N increases, does the convergent sequence of interpolating polynomials corresponding to the 
optimal solutions of Problem |2] converge to an optimal solution of Problem [T]? Namely, 

lim {Ijsix\Inu^) — ix*,u*) 

TV— s-oo 

Remark 2: It is possible that Problem [T] has more than one optimal solution, i.e., there is more than one solution 
with the same optimal cost J{x*,u*) = J*. Therefore, to show that the sequence of discrete solutions converges 
to an optimal solution, we can instead show that the cost of the discrete solution, J, converges to the optimal cost 
J*. 

Previous work has been done in the area of convergence of the pseudospectral method and we aim to augment this 
literature with several key insights that make convergence results appUcable to a wider class of systems and relax the 
conditions on which the current proofs are based. Rather complete analysis has been done for the class of nonlinear 
systems which can be feedback linearized, including convergence rates ll36l . We show below that ensemble quantum 
systems of interest do not fall within the class of feedback linearizable systems. Work has also included general 
nonlinear systems, but with the assumption that the solutions of the algebraic nonlinear programming problem have 
a limit point (i.e., have a convergent subsequence) ll37l . In the language used above, this is very close to assuming 
"Convergence", which in this presentation we relax and prove Feasibility, Convergence, and Consistency directly. 
Finally, we examine the convergence of the multidimensional pseudospectral method as applied to ensemble optimal 
control problems. In what follows we consider first the convergence of the standard pseudospectral method and 
then discuss the convergence of the ensemble case. 

We first observe that ensemble control systems of interest are not feedback linearizable ll38l . which motivates 
a need for a more general convergence proof. Consider the bilinear Bloch equations in (fT9] l without variation in 
rf inhomogeneity (i.e., e — 1). The ability to feedback linearize a general nonlinear system is given by the Lie 
algebra generated by the drift and control vector fields (the conditions on this algebra must hold for each control 
term individually; here we consider the case for u). In particular, the terms ad° = fly, ad^J^^^ ily — -^uiflx, 
ad^j2^57j, = —Lj'^fly, . . . , and, 

where k = 1,2,..., and uj is any value in the interval C M. It is clear that this Lie algebra, with increasing 
powers of the parameter lu, is never closed. Therefore, the span of the appropriate Lie brackets is not involutive, 
which indicates that such a system is not feedback linearizable. 

A. Empirical Convergence 

The orthogonal polynomials of the pseudospectral method provide spectral convergence rates similar to Fourier 
series approximations for periodic functions, which can easily be seen in practice. Figure |4] shows the rapid 
convergence of the method in both the discretization (time) and sampling (parameter) dimensions for a broadband 
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Fig. 4. The characteristic rapid convergence of the multidimensional pseudospectral method for a tt/2 pulse designed to perform the state 
transfer Af(0,aj, 1) = (0 1)' -> M{T,uj, 1) = (1 0)', with _B = 1 and T = 1. Average terminal values of Mx{T,u}, 1) are shown for 
various choices of A' and Nuj. 

7r/2 pulse maximizing the terminal x value across the ensemble. As the order of discretization (N) and/or sampling 
{Ns) increase, the method yields an objective {f{T) ~ Mx{T, u, 1)) that converges to the maximum value of unity. 
The low order of approximation is a characteristic of the orthogonal approximations at the heart of the numerical 
method. Although such empirical figures are convincing, we now show this convergence in a more rigorous fashion. 

B. Theoretical Preliminaries 

The results in this section will provide the foundation on which we can analyze the feasibility, convergence, 
and consistency of the pseudospectral approximation method for optimal control problems. We begin by presenting 
several key established results in polynomial approximation theory and the natural vector extensions. With these 
identities, we are able to then prove feasibility and convergence. We define an optimal solution to Problem [T] as 
any feasible solution that achieves the optimal cost J{x*,u*) = J*. We use this definition of an optimal solution 
within the subsequent preliminaries and the main result. 

Remarks.- Given Problem [T] x E H"{Q). Since x{t) exists and / G C""\ all the derivatives a;^'') G C°, 
V A; = 0, 1, . . . , a exist and are square integrable on the compact domain fl, x^*"'' G L'^{Q). Therefore, x G H"{Q). 

Lemma 1 (Interpolation Error Bounds 4251/ . p. 289): If /i G H°'{fl), the following hold with ci, C2, C3, c > 0. 

(a) The interpolation error is bounded, 

\\h-lNh\\2<ciN-''\\h\\^^y 

(b) The error between the exact derivative and the derivative of the interpolation is bounded, 

\\h-VNh\\2<C2N^-''\\h\\^c)- 

The same bound holds for the discrete norm, 

\\h-VNh\\N <C3N^-"\\h\\i^o.)- 
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(c) The error due to quadrature integration is bounded, 

N 



„1 iV 

/ h{t)dt - h{tk)wk 
•^-1 fc=o 



< 



c7v-"||;.||(„), 



where tk is the fc* LGL node and Wk is the corresponding fc* weight for LGL quadrature as defined in (ITZt . 

Lemma 2: If ft, e H^{il), i.e., an n-vector valued Sobolev space, h — {hi h2 ... hn)', hi g H°'{n), i = 
l,2,...,n. 

(a) The vector-valued extension of Lemma flftl is. by the triangle inequality on the Lf^{n) norm, 

n n 
\\h-lNhh < E \\h^-INh^\\2 <Y.C,N-"\\h,\\^^y 
i=l i=l 

(b) Similarly, [30 can be extended, 

n n 
1=1 i=l 

which again also holds for the discrete -L^j(ri) norm. 

Proposition 1 (Feasibility): Given a solution {x, u) of ProblemlT] then Problem|2]has a feasible solution, (x, w), 
which are the corresponding interpolation coefficients. 

Proof: Given the feasible solution (a;, u), let {Inx,Inu) be the polynomial interpolation of this solution at the 
LGL nodes. Our aim is to show that the coefficients of this interpolation satisfy (|22]|-(|24]| of Problem |2] Consider 
the constraints imposed by the dynamics in ( |22] |. Because the discrete norm is evaluated only at the interpolation 
points, 

\\f{lNX,lNu) - Vnx\\n = \\f{x,u) - Vnx\\n = \\x - r>Ara;||Ar < CdiV^"" 

where the last step is given by Lemma QHl Therefore, the interpolation coefficients {x,u) satisfy the dynamics of 
Problem|2]in (l22l i. We can easily show that the path constraints are also satisfied because g{x{t), u{t)) < for all 
t e by (Ull. Since this holds for all t £ fl, it also holds for all LGL nodes tk G T^'^^, i.e., 

g{xk,Uk) = g{x{tk),u{tk)) < 0, 

which gives (|24] |. The endpoint constraints are trivially satisfied by the definition of interpolation and the presence 
of interpolation nodes at both endpoints. Therefore, {x,u) is a feasible solution to Problem |2] ■ 

Proposition 2 (Convergence): Given the sequence of solutions to Problem |2l {{x,u)}n, then the sequence of 
corresponding interpolation polynomials, {{Inx,Inu)}, has a convergent subsequence, such that 

lim (Ijv, a; , Iat w) = (loo a; , loo , 

Nj — f oo 

which is a feasible solution to Problem [T] 

Proof: Given that {x, u) is a feasible solution of Problem |2] it satisfies ([22]l-(l24ll. Our goal is to show (i) 
that the sequence of solutions, {{Inx,Inu)}n, has a convergent subsequence and (ii) that the limit point of this 
function subsequence is a feasible solution of Problem [T] satisfying ©-(ISll. 
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(i) The sequence {I^x}, is a sequence of polynomials on a compact domain, therefore, J^x G H"{il). With 
the boundedness of the interpolating polynomials and the compactness of O, Rellich's Theorem (see Appendix |B]i 
states there is a subsequence {Ip^.x} which converges in H"^^{il). The same is true for the control interpolating 
polynomial. Therefore, there exists at least one limit point of the function sequence {{Inx,Xnu)} which we denote 

(ii) Explicitly writing out the calculation of the discrete norm in (l22l i gives 

1/2 



' N n \ ^/^ 



In the limit, because / is continuous, 

lim (fi{lNX,lNu) - VNx){tk) = (fi{IocX,IocU) - {Iooxy){tk) = 0, 

therefore, 

at 

which states that [XooX^Ioau) satisfies the dynamics in ^ at the interpolation nodes. Moreover, as ^> oo, the 
LGL nodes G r'"'^^ are dense in Vt, which further shows that {IaoX^Icx,u) satisfies the dynamics of Problem 
[T] at all points on the interval Vl. Similarly, one can prove that this solution satisfies the path constraints because 
the LGL nodes become dense in O as ^ cxa and g{xk,Uk) = g{x{tk), u{tk)) < at all LGL nodes. Again, the 
endpoint constraints are met exactly because the LGL grid has nodes at the endpoints. ■ 

Lemma 3: Given {x, u), where x e H^{il), u G H^^{n), and the corresponding interpolation coefficients, {x, u), 
then the error in the continuous and discrete cost functionals defined in (O and d^Tt . respectively, due to interpolation 
is given by, 

|J(x,w) - J(x,u)| < cN-". 

Remark 4: Notice that {x, u) and (x, u) are not required to be a feasible solutions to Problem [Hand |2] respectively. 
This result characterizes the error due to interpolation. 

Proof: From (|6]l and (l22T i since (p{x{l)) = (p{xn), 

N 

\J{x, u) — J(x, u)\ = 

' ' ' fe=0 

Since C G C" with respect to both the state and control, x G if," (57) and u G H^{il), the composite function 
C{t) = C{x{t),u{t)) G iJ"(51). Let Ck = C{xk,Uk)- Substituting these definitions and employing Lemma \W\ we 
obtain 



»1 iV 

/ C{x, u)dt - C{xk,Uk)wk 



.1 N 

/ c{t)dt 

— 1 1 n 



kWk 

k=0 



<ciV-"||£(t)||(„). 



Since jC g iJ"(ri), is bounded and the result follows. 
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V. Main Result 

Theorem 1 (Consistency): Suppose Problem[T]has an optimal solution (x*, u*). Given a sequence of optimal solu- 
tions to Problem|2] {{x\ ut)}jv, then the corresponding sequence of interpolating polynomials, {{Imx\Xmu^)}n, 
has a limit point, {Xoox\looU^) which is an optimal solution to the original optimal control problem. 

Proof: We break the proof into four sections, employing the results from the previous section. 

(i) By Proposition [T] since is a solution to Problem [T] then for each choice of N, the corresponding 
interpolation coefficients, [x* ,u*), are a feasible solution to Problem |2] By the definition of optimality of {x\u^), 

J{x\v)) <J{x*,u*). (26) 

( ii) By Proposition |2] the limit point of the polynomial interpolation of the discrete optimal solution to Problem |2] 
limiv^oo(^Afa;^,Ijvw^) = (looa^^j^TooW^), is a feasible solution of Problem[T| Therefore, we have, by the definition 
of the optimality of {x*,u*) and the continuity of J, 

J{x*,u*)< lim JilNX^lNu') = J{Xocx\l.^u'). (27) 

N^oo 

(Hi) Using Lemma |3] we can bound the error in the cost between the optimal solution of Problem [T] {x* ,u*), and 
the corresponding interpolating coefficients, {x*,u*), as 

\J{x*,u*) - J(x*,r)| < ciN-°'. (28) 

Similarly, we can bound the error in the cost between the optimal solution of Problem|2] (x^, u^), and the polynomial 
interpolation of this solution, (Ijva;^,lArM^), as 

\J{Inx^,Inu^) - J{x\u'')\ < C2N-". (29) 

Recall that Lemma [3] does not require (I.mx\Xmu^) to be a feasible solution of Problem [T] From (l28l l and (|29] l, 

lim Jix* ,u*) = J{x* ,u*), (30) 
lim \J{Xnx\Xnu^) ~ J{x\v))\ ^Q. (31) 

AT— ^00 

(iv) We are now ready to assemble the various pieces of this proof. Combining ( l30l l and ( l26l l we have, 

lim J{x\v))< lim J{x\u*) ^ J{x\u*). 

AT— >oo TV— s-oo 

Adding the result from (l27t . 

lim J{x\v)) < J{x*,u*) < lim J{Xnx\Xnu^). (32) 

N-^oo N—i-oo 

Since the difference between the left and right sides, as given by (l3Tl i. decreases to zero as — 00, the quantities 

J(a;t,ut) and J{Xis!x\Xmu^) converge to J{x*,u*). In particular, 

0< lim \j(x*,u*)-J(x\u'')] < lim \J(Xnx\Xnu'')-J{x\u'')]=0. 
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Thus the optimal discrete cost J{x\ u^) of Problem |2] and the continuous cost J{2nx'^ jJnu^) of the corresponding 
interpolation polynomials converge to the optimal cost J{x*,u*) of Problem [T] Moreover, (looa^^j^^oo^*^) is a 
feasible solution to Problem [T] and achieves the optimal cost. Therefore, (Iooa;^,XooM^) is an optimal solution to 
Problem m ■ 
Remark 5 (Ensemble Extension): The nature in which the ensemble extension enters into the multidimensional 
pseudospectral method makes it straightforward to extend this convergence proof to the ensemble case. Section 
III-DI showed the simplicity of the derivative term in multidimensional sampling with equation ( fTSl ). Similarly, in the 
ensemble case, the constraints corresponding to the dynamics (l22l) operate entirely in parallel for different parameter 
values. The additional integration in the cost function over the parameter domain, as in (|4) adds another layer of 
quadrature approximation that can be shown to converge with arguments similar to those presented above. 

VI. Conclusion 

In this work we have presented a cohesive perspective and methodology for optimal control of inhomogeneous 
ensembles, as particularly motivated by compelling problems in quantum control and extendable to both parame- 
terized systems in, for example, neuroscience El and uncertain systems throughout science and engineering. Such 
systems are mathematically characterized by considering a parameterized family of differential equations indexed 
by a parameter vector that shows variation. Applying this rigorous framework prompts us to solve the corresponding 
optimal control problems with computational methods of particular form. The notion of polynomial approximation 
entering into the controllability analysis of the Bloch equations indicates that a modified pseudospectral method 
is a prime candidate. The method has natural extensions which we develop to model ensemble variation. This 
direct collocation method transforms the continuous-time optimal control problem into an algebraic nonlinear 
programming problem, which we show to be effective in a variety of applications. We supplied additional and more 
general arguments for the convergence of this method, in particular relaxing several assumptions and discussing the 
convergence characteristic of the multidimensional pseudospectral method for optimal ensemble control. 



Appendix A 
The Dimensionless Bloch Equations 

The Bloch equations without relaxation, Ai ~ A4 x jBf^ff, utilizes the classical description of interacting 
electromagnetic forces, where A4 is the spin magnetization vector, 7 is the gyromagnetic ratio, the effective 
externally applied field is Z?eff = {Bi cos{ujot + (p), Bi sm{ujot + </>), B^)', Bi{t) and So are the amplitudes of the 
applied fields in the transverse plane and z direction respectively, and is the phase angle H]. Conventionally, 
the fields are given as frequencies ^B^fs — (wi^, wij^, wo) and measured in units of Hertz. Using the generators of 
rotation. 
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the Bloch equations are be given by 

d 



dt 



M{t) ^ LOq^z + ^^ly{t)^y + i^lx{t)^x M{t). 



(33) 



If we consider variation in the apphed electromagnetic fields Bq and Bi, we can express (|33T l in matrix form. 



d_ 

di 



Mxit,uj,e) 
My{t,uj,e) 
Mz{t,oj,e) 



-(wo+w) 
Wo + w 
—eBi sm{LUot + (jj) eBi cos{uj{)t + 



eBi sin(a;ot + 0) 
—eBi cos{ujot + (j)) 




My{t,uj,e) 



where lu G [— /3,/3] and e G [1 — S,l + S], < 6 < 1. For calculation and computation, it is useful to transform 
the Bloch equations into the so-called rotating frame and normalize the system by a nominal pulse amplitude A 
to yield a dimensionless equation. Solutions based on the dimensionless equation can then be scaled for a specific 
choice of nominal amplitude. Consider a transformation M ~ exp{~ujoilzt)A4. In addition we scale time with 
T = At. It is straightforward to show that the new state equation is given by. 



dr 

where t e [0, AT x 2tt], oj e [-B, B], B = /3/A, and 

iB,{t/A) 



u{t) = 



■ cos 



^(r)^^^^4^sin(0(.M)), 



A V-w V V / ^ 

(aU dimensionless). Note the 2tt factor in the time scaling is introduced to convert from units of Hertz to radi- 
ans/second. Designing the time-varying controls u{t) and v{t) is equivalent to the original design of amplitude 
Bi{t) and phase 4'{t). 

Appendix B 
Rellich's Theorem 

Theorem 2 (Rellich's Theorem [39'], p. 272): Suppose {/a,} is a sequence in such that 

(i) supfc 1 1 /fell (a) < oo, and 

(ii) the /fe's are all supported in a fixed compact set V. 

Then there is a subsequence {fk } which converges in for all f3 < a. 
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